/*************************************************************/
/*   Source for a test of a new elasticity tensor in MOOSE   */
/*                  Status: not working                      */
/*                                                           */
/* Copyright (C) 2017, 2018  eDgar                           */
/*                                                           */
/* This program is free software: you can redistribute it    */
/* and/or modify it under the terms of the GNU General       */
/* Public License as published by the Free Software          */
/* Foundation, either version 3 of the License, or (at your  */
/* option) any later version.                                */
/*                                                           */
/* This program is distributed in the hope that it will be   */
/* useful, but WITHOUT ANY WARRANTY; without even the        */
/* implied warranty of MERCHANTABILITY or FITNESS FOR A      */
/* PARTICULAR PURPOSE.  See the GNU General Public License   */
/* for more details.                                         */
/*                                                           */
/* You should have received a copy of the GNU General Public */
/* License along with this program.  If not, see             */
/* <http://www.gnu.org/licenses/>.                           */
/*************************************************************/

// We need to get the specification of how to call the
// functions which will be defined below. That is the main
// point of the header (.h) file
#include "spider.h"

// Populate the template which was sketched in the header
// file. The parameters given here can be used by the user
// in the input file.
template <>
InputParameters
validParams<spider>()
{
  InputParameters params = validParams<ComputeElasticityTensorBase>();

  // Description which is shown in the auto-generated
  // documentation
  params.addClassDescription("DO NOT USE!!\n"
                             "spider defines an "
                             "elasticity tensor for the "
                             "St. Venant-Kirchhoff material. "
                             "It is meant as a guide on the "
                             "implementation of (elasticity "
                             "tensors of) new materials.\n"
                             "\\("
                             "\\Psi(\\tensor{\varepsilon})"
                             "= \\frac{1}{2}\\,"
                             " \\lambda\\,"
                             "  (\\Tr{"
                             "    \\tensor{\varepsilon}"
                             "  })^{2} +"
                             "\\mu\\,"
                             " \\tensor{\varepsilon}:"
                             "  \\tensor{\varepsilon}"
                             "\\)"
                             );

  // State that _lambda is a required parameter
  params.addRequiredParam<Real>("lambda",
                        "Lame's first constant for the material.");
  // Give a default value to \mu (set -1; if the value is
  // not set to a positive value by the user, there is a
  // check below which throws an error)
  params.addParam<Real>("mu",
                        -1.0,
                        "Lamé second constant for the material");

  return params;
}

// Notice that the parameters become constants
// when we make this constructor
spider::spider(const InputParameters & parameters)
  : ComputeElasticityTensorBase(parameters),

    // When this function is called, it is because the
    // application is compiled and executed with an input
    // file (.i).
    //
    // First, make sure that the parameters were given in
    // the .i (input) file. Note that we called it "lambda"
    // in the template above, but the variable in the .h
    // file is _lambda. For the user filling the .i (input)
    // file, it will be `lambda'. For calculations in the
    // code, it will be `_lambda'.
    _lambda_set(parameters.isParamValid("lambda")),
    _mu_set(parameters.isParamValid("mu")),

    // Get parameters from the input file. The following
    // line is an implicit conditional. It checks if the
    // Boolean value is True, in which case, it assigns the
    // value after the `?' (and before `:'). If the value is
    // False, then the variable gets the value indicated
    // after `:'. In this case, if `lambda' was set in the
    // .i (input) file, it gets the value from
    // there. Otherwise, it sets it to -1.
    _lambda(_lambda_set ? getParam<Real>("lambda") : -1),
    // The same for `mu'.
    // _mu(_mu_set ? getMaterialProperty<Real>("mu") : -1)
    _mu(_mu_set ? getParam<Real>("mu") : -1)
{
  // St. Venant-Kirchhoff material model (example; not
  // useful for any practical application).

  // Now we are going to initialise the tensor

  // Make sure that the values are not both negative or zero
  // Holzapfel, G. A., Nonlinear solid mechanics : a
  // continuum approach for engineering (2000), : Wiley. p.250
  if (_lambda < 0.0 || _mu < 0.0)
    mooseError("The Lamé constants must be positive in material '"
               + name() + "'.");

  // The St. Venant-Kirchhoff model is defined as
  //
  // LaTeX:
  // \[{\cal C}_{ijkl} = \lambda\,\delta_{ij}\,\delta_{kl}
  //                     + \mu\,\left(
  //                         \delta_{ik}\,\delta_{jl}
  //                         + \delta_{il}\,\delta_{jk}
  //                       \right) \]
  //
  // You have to make sure that whatever operatations
  // (including addition, subtraction, etc. work as you
  // think). In this case, check that the definition of
  // RankFourTensor has the required operations.
  // [[moose/framework/src/utils/RankFourTensor.C]]

  // // Identity tensor (top left 3x3 diagonal of a fourth order
  // // tensor in Voigt's notation; not used)
  // // 1  0  0  0  0  0
  // // 0  1  0  0  0  0
  // // 0  0  1  0  0  0
  // // 0  0  0  0  0  0
  // // 0  0  0  0  0  0
  // // 0  0  0  0  0  0
  // // initIdentity() in
  // // [[moose/framework/src/utils/RankFourTensor.C]]
  // const RankFourTensor tensI = RankFourTensor::initIdentity;

  // LaTeX:
  // \( \left(\delta_{ik}\,\delta_{jl}
  //    + \delta_{il}\,\delta_{jk}\right) \)
  //
  // 1.0  0.0  0.0  0.0  0.0  0.0
  // 0.0  1.0  0.0  0.0  0.0  0.0
  // 0.0  0.0  1.0  0.0  0.0  0.0
  // 0.0  0.0  0.0  0.5  0.0  0.0
  // 0.0  0.0  0.0  0.0  0.5  0.0
  // 0.0  0.0  0.0  0.0  0.0  0.5
  // initIdentitySymmetricFour() in
  // [[moose/framework/src/utils/RankFourTensor.C]]
  //
  // initIdentitySymmetricFour
  // gives (de_ij de_kl + de_il de_jk)/2 -- andrew.wilkins
  //
  // It's possible that someone might use this with crazy
  // non-symmetric stresses and strains and
  // initIdentitySymmetricFour essentially symmetrises the
  // strain before forming the stress. -- andrew.wilkins
  //
  // See also [BelytschkoLiuMoranElkhodary14, p.240]
  const RankFourTensor dik_djl = RankFourTensor::initIdentitySymmetricFour;

  // Identity matrix
  // 1  0  0
  // 0  1  0
  // 0  0  1
  // initIdentity() in
  // [[moose/framework/src/utils/RankTwoTensor.C]]
  // const RankTwoTensor matI = RankTwoTensor(RankTwoTensor::initIdentity);
  const RankTwoTensor dij(RankTwoTensor::initIdentity);

  // LaTeX:
  // \( \delta_{ij} \, \delta_{kl} \)
  //
  // 1.0  1.0  1.0  0.0  0.0  0.0
  // 1.0  1.0  1.0  0.0  0.0  0.0
  // 1.0  1.0  1.0  0.0  0.0  0.0
  // 0.0  0.0  0.0  0.0  0.0  0.0
  // 0.0  0.0  0.0  0.0  0.0  0.0
  // 0.0  0.0  0.0  0.0  0.0  0.0
  const RankFourTensor dij_dkl = dij.outerProduct(dij);

  // // Another way for dij_dkl
  // // To make sure that nothing is out of place, fill the
  // // whole tensor with zeros
  // // [[moose/framework/include/utils/RankFourTensor.h]]
  // dij_dkl.zero();
  // // Set the corresponding values
  // // [[moose/framework/src/utils/RankFourTensor.C]]
  // dij_dkl.fillGeneralIsotropic(1, 0, 0);

  // This is the right definition of St. Venant-Kirchhoff
  // Knowing that the * operator is defined as we expect
  // between a tensor and a scalar, as well as the +
  // operator between the two RankFourTensor
  // [[moose/framework/src/utils/RankFourTensor.h]]
  // , from the way the St. Venant-Kirchhoff material is
  // defined:
  _Cijkl = _lambda * dij_dkl + 2 * _mu * dik_djl;

  // Using the symmetric isotropic function provided by
  // MOOSE to produce the St. Venant-Kirchhoff tensor:
  //
  // λ + 2 µ     µ        µ     0.0000  0.0000  0.0000
  //    µ     λ + 2 µ     µ     0.0000  0.0000  0.0000
  //    µ        µ     λ + 2 µ  0.0000  0.0000  0.0000
  // 0.0000   0.0000   0.0000      µ    0.0000  0.0000
  // 0.0000   0.0000   0.0000   0.0000     µ    0.0000
  // 0.0000   0.0000   0.0000   0.0000  0.0000     µ
  // Thanks to Stephanie Pitts
  // [[http://mooseframework.org/docs/moose_docs/site/
  //   documentation/systems/Materials/tensor_mechanics/
  //   ComputeIsotropicElasticityTensor/index.html]]
  // See RankFourTensor.C::fillSymmetricIsotropic,
  //     RankFourTensor::fillSymmetricIsotropicEandNu
  //
  // (clear the values first; just to make sure and show
  // this function)
  _Cijkl.zero();
  // [[moose/framework/src/utils/RankFourTensor.C]]
  _Cijkl.fillSymmetricIsotropic(_lambda, _mu);

  // // Another way to do it is by setting the values
  // // directly (note that the following is incomplete)
  // // Thanks to Daniel Schwen:
  // //    _Cijkl(0,0,0,0). Rank Four needs 4 indices (and
  // //    for that you need the round bracket operator)!
  // for (uint i = 0; i < 3; i++)
  //   _Cijkl(i, i, i, i) = 1.0;
  //
  // for (uint i = 0; i < 3; i++)
  //   {
  //     for (uint j = 0; j < 3; j++)
  //       {
  //         for (uint k = 0; k < 3; k++)
  //           {
  //             for (uint l = 0; l < 3; l++)
  //               {
  //                 _Cijkl(i, j, k, l) = (i==j)&&(l==k);
  //               }
  //           }
  //       }
  //   }

  // All tensors created by this class are always isotropic
  // (ComputeFiniteStrainElasticStress can only be used with
  // elasticity tensors that guarantee isotropy).
  if (_Cijkl.isSymmetric())
    issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);

  // If showing the values of the tensor is required, this
  // is a hackish solution
  FILE *fp;

  fp = fopen("test_spider.txt", "w");

  /*
  \begin{bmatrix}
  C_{1111} & C_{1122} & C_{1133} & C_{1123} & C_{1131} & C_{1112} \\
  C_{2211} & C_{2222} & C_{2233} & C_{2223} & C_{2231} & C_{2212} \\
  C_{3311} & C_{3322} & C_{3333} & C_{3323} & C_{3331} & C_{3312} \\
  C_{2311} & C_{2322} & C_{2333} & C_{2323} & C_{2331} & C_{2312} \\
  C_{3111} & C_{3122} & C_{3133} & C_{3123} & C_{3131} & C_{3112} \\
  C_{1211} & C_{1222} & C_{1233} & C_{1223} & C_{1231} & C_{1212}
  \end{bmatrix}
  */

  // % 10.4e : print a number with exponential notation (e),
  // leaving a blank space if it's non-negative (space
  // after %), with 3 decimals (.3) and reserving 10
  // characters for the number with its mantissa.
  // Print the elasticity matrix in Voigt's notation.
  fprintf(fp,
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "% 10.4e % 10.4e % 10.4e % 10.4e % 10.4e % 10.4e \n"
          "\n",
          _Cijkl(0,0,0,0), _Cijkl(0,0,1,1), _Cijkl(0,0,2,2), _Cijkl(0,0,1,2), _Cijkl(0,0,2,0), _Cijkl(0,0,0,1),
          _Cijkl(1,1,0,0), _Cijkl(1,1,1,1), _Cijkl(1,1,2,2), _Cijkl(1,1,1,2), _Cijkl(1,1,2,0), _Cijkl(1,1,0,1),
          _Cijkl(2,2,0,0), _Cijkl(2,2,1,1), _Cijkl(2,2,2,2), _Cijkl(2,2,1,2), _Cijkl(2,2,2,0), _Cijkl(2,2,0,1),
          _Cijkl(1,2,0,0), _Cijkl(1,2,1,1), _Cijkl(1,2,2,2), _Cijkl(1,2,1,2), _Cijkl(1,2,2,0), _Cijkl(1,2,0,1),
          _Cijkl(2,0,0,0), _Cijkl(2,0,1,1), _Cijkl(2,0,2,2), _Cijkl(2,0,1,2), _Cijkl(2,0,2,0), _Cijkl(2,0,0,1),
          _Cijkl(0,1,0,0), _Cijkl(0,1,1,1), _Cijkl(0,1,2,2), _Cijkl(0,1,1,2), _Cijkl(0,1,2,0), _Cijkl(0,1,0,1)
          );

  // You can print all the components like this
  for (int i = 0; i < 3; i++)
    {
      for (int j = 0; j < 3; j++)
        {
          for (int k = 0; k < 3; k++)
            {
              for (int l = 0; l < 3; l++)
                {
                  fprintf(fp, "%d%d%d%d: % 10.4e\n", i + 1, j + 1, k + 1, l + 1, _Cijkl(i, j, k, l));
                }
            }
        }
    }

  // Remember to *ALWAYS* close the file pointer (descriptor)
  fclose(fp);

  // You can get a second order tensor from two vectors like this:
  // 1. Declare the vector this (either works)
  // const TypeVector<Real>;
  const RealVectorValue vec;     // RealVectorValue is
                                // declared in
                                // vector_value.h (libMesh)
  // 2. Declare the second order tensor
  RankTwoTensor A;              // This cannot be a constant
  // 3. Do the outer product
  A.vectorOuterProduct(vec, vec);
  // 4. You can also get the outer product of two second order tensors
  const RankFourTensor tens4 = A.outerProduct(A);
}

void
spider::computeQpElasticityTensor()
{
  // This assigns an elasticity tensor at each quadrature
  // point (_qp) to each of the elements of the material.
  _elasticity_tensor[_qp] = _Cijkl;
}
